#!/usr/bin/perl -w
# Given a specified sequence and a specified set of sequences, first finds the
# percent identity between the single sequence and each member of the sequence
# set, then gets the mean and standard deviation of these percent IDs.
use strict;
use warnings;

if (@ARGV != 2) {
  print "USAGE: $0 protein.fasta output.fasta\n";
  exit(1);
}

# assumes sequences are all of equal length
my ($fastaSeedSeqFile, $fastaFile) = @ARGV; 

my %codonTable = (
 "AAA" => "K", "AAC" => "N", "AAG" => "K", "AAT" => "N",
 "ACA" => "T", "ACC" => "T", "ACG" => "T", "ACT" => "T",
 "AGA" => "R", "AGC" => "S", "AGG" => "R", "AGT" => "S",
 "ATA" => "I", "ATC" => "I", "ATG" => "M", "ATT" => "I",

 "CAA" => "Q", "CAC" => "H", "CAG" => "Q", "CAT" => "H",
 "CCA" => "P", "CCC" => "P", "CCG" => "P", "CCT" => "P",
 "CGA" => "R", "CGC" => "R", "CGG" => "R", "CGT" => "R",
 "CTA" => "L", "CTC" => "L", "CTG" => "L", "CTT" => "L",

 "GAA" => "E", "GAC" => "D", "GAG" => "E", "GAT" => "D",
 "GCA" => "A", "GCC" => "A", "GCG" => "A", "GCT" => "A",
 "GGA" => "G", "GGC" => "G", "GGG" => "G", "GGT" => "G",
 "GTA" => "V", "GTC" => "V", "GTG" => "V", "GTT" => "V",

 "TAA" => "X", "TAC" => "Y", "TAG" => "X", "TAT" => "Y",
 "TCA" => "S", "TCC" => "S", "TCG" => "S", "TCT" => "S",
 "TGA" => "X", "TGC" => "C", "TGG" => "W", "TGT" => "C",
 "TTA" => "L", "TTC" => "F", "TTG" => "L", "TTT" => "F" 
);

# expecting single fasta protein
my %seedSeq = getFastaFrom ($fastaSeedSeqFile);
if ((values %seedSeq) == 0) {
  print "Need at least one sequence in $fastaSeedSeqFile.\n";
  exit(1);
}

my $seed = (values (%seedSeq))[0];
$seed = uc $seed;

# expecting multiple fasta dna sequences (hence using translate procedure)
my %progenySeqs = getFastaFrom ($fastaFile);

if ((values %progenySeqs) == 0) {
  print "Need at least one sequence in $fastaFile.\n";
  exit(1);
}

my %generatedSeqs = translate (\%codonTable, %progenySeqs);
my @fracIDs;

# get the identities between each sequence and the seed sequence, then put
# them in @fracIDs
foreach my $progeny (values %generatedSeqs) {
  $progeny = uc $progeny;
  my $identity = fractionID ($seed, $progeny);
  push @fracIDs, $identity;
}

my $mean = mean (@fracIDs);
my $stddev = stdDeviation ($mean, @fracIDs);

printf "Mean Percent ID +/- Standard Deviation: %.1f +/- %.1f\n",100*$mean,100*$stddev;

# given a mean and a data sample, returns the standard deviation
sub stdDeviation {
  my ($mean, @observed) = @_;
  my $sum=0;
  foreach my $datum (@observed) {
    $sum += ($datum - $mean)**2;
  }
  my $variance = $sum/@observed;
  return sqrt ($variance);
}

# given a data sample, returns the mean
sub mean {
  my (@observed) = @_;
  my $sum=0;
  foreach my $datum (@observed) {
    $sum += $datum;
  }
  return $sum/@observed;
}

# returns fraction ID of a pair of sequences (ignoring gaps, and assuming 
# strings have equal length)
sub fractionID {
  my ($str1,$str2) = @_;
  $str1 =~ s/\s//g;
  $str2 =~ s/\s//g;
  # split into individual characters
  my @seq1 = split (//, $str1);
  my @seq2 = split (//, $str2);
  if (length(@seq1) != length(@seq2)) {
    die "Error: found sequences of different lengths. Exiting...";
  }

  my $matchCount=0;
  for (my $i=0; $i<@seq1; $i++) {
    if ($seq1[$i] eq $seq2[$i]) {
      $matchCount++;
    }
  }

  return $matchCount/@seq1;
}

sub translate {
  my ($codonTable, %dnaFasta) = @_;
  my %codonTable = %{$codonTable};
  my %proteinFasta;
  foreach my $label (keys %dnaFasta) {
    my $seq = $dnaFasta{$label};
    $seq =~ s/^.*?([ACGT]{3}.*?)[ACGT]{3}.*$/$1/;
    $seq = uc $seq; 
    my @codons = ( $seq =~ m/.{3}/g );
    my $protein = "";
    for (my $i=0; $i<@codons; $i++) {
      $protein .= $codonTable{$codons[$i]};
    }
    $proteinFasta{$label} = $protein; 
  }
  return %proteinFasta;
}

# including sequence labels
sub getFastaFrom {
  my ($fastafile) = @_;
  my %fasta;
  open (INPUT, $fastafile) or die "Can't open $fastafile: $!";

  my $currentSequence;
  while( my $line = <INPUT> ) {
    chomp $line;
    if( $line =~ /^\s*$/ ) {
      next;
    } elsif( $line =~ /^>/ ) {
      $currentSequence = $line;
      $fasta{$currentSequence} = "";  
    } else {
      # concatenate the line to the current sequence if it is not a label
      $line =~ s/\s//g;
      $fasta{$currentSequence} .= $line;
    }
  }
  return %fasta;
}

